%% it computes transition matrices for the HJB equation

Db_f = zeros(Ntot,Ntot);
for b=1:N_b-1
    for z=1:N_z
        Db_f((b-1)*N_z+z,(b-1)*N_z+z) = -1/db;
        Db_f((b-1)*N_z+z,b*N_z+z)     =  1/db;
    end
end
Db_f=sparse(Db_f);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Db_b = zeros(Ntot,Ntot);
for b=1:N_b-1
    for z=1:N_z
        Db_b(b*N_z+z,b*N_z+z)     =  1/db;
        Db_b(b*N_z+z,(b-1)*N_z+z) = -1/db;
    end
end
 Db_b=sparse(Db_b);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tz_f = zeros(Ntot,Ntot);
for b=1:N_b
    for z=1:N_z-1
        Tz_f((b-1)*N_z+z,(b-1)*N_z+z)    = -1/dz;
        Tz_f((b-1)*N_z+z,(b-1)*N_z+z+1)  =  1/dz;
    end
end
Tz_f=sparse(Tz_f);

Tz_b = zeros(Ntot,Ntot);
for b=1:N_b
    for z=2:N_z
        Tz_b((b-1)*N_z+z,(b-1)*N_z+z)    =  1/dz;
        Tz_b((b-1)*N_z+z,(b-1)*N_z+z-1)  = -1/dz;
    end
end
Tz_b=sparse(Tz_b);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tzz  = zeros(Ntot,Ntot);%this is assuming there is similar mass around z=zmin and z=zmax
for b=1:N_b
    Tzz((b-1)*N_z+1,(b-1)*N_z+1)        = -1/dz^2;
    Tzz((b-1)*N_z+1,(b-1)*N_z+2)        =  1/dz^2;
    for z=2:N_z-1
        Tzz((b-1)*N_z+z,(b-1)*N_z+z-1)  =  1/dz^2;
        Tzz((b-1)*N_z+z,(b-1)*N_z+z)    = -2/dz^2;
        Tzz((b-1)*N_z+z,(b-1)*N_z+z+1)  =  1/dz^2;
    end
    Tzz((b-1)*N_z+N_z,(b-1)*N_z+N_z-1)  =  1/dz^2;
    Tzz((b-1)*N_z+N_z,(b-1)*N_z+N_z)    = -1/dz^2;
end
Tzz=sparse(Tzz);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dbb  = zeros(Ntot,Ntot);%this is assuming there is no mass aaround b=bmin and b=bmax
for z=1:N_z
    for b=2:N_b-1
        Dbb((b-1)*N_z+z,(b-2)*N_z+z)  =  1/db^2;
        Dbb((b-1)*N_z+z,(b-1)*N_z+z)  = -2/db^2;
        Dbb((b-1)*N_z+z,    b*N_z+z)  =  1/db^2;
    end
end
Dbb=sparse(Dbb);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dbz  = zeros(Ntot,Ntot);
for z=2:N_z-1
    for b=2:N_b-1
        Dbz((b-1)*N_z+z,(b  )*N_z+z+1)=  1/(4*db*dz);
        Dbz((b-1)*N_z+z,(b-2)*N_z+z+1)= -1/(4*db*dz);
        Dbz((b-1)*N_z+z,    b*N_z+z-1)= -1/(4*db*dz);
        Dbz((b-1)*N_z+z,(b-2)*N_z+z-1)=  1/(4*db*dz);
    end
end
Dbz=sparse(Dbz);